Grizzly bears (Ursus arctos) in the Greater Yellowstone Ecosystem of the US Rocky Mountains have recently increased in numbers, but remain vulnerable due to isolation from other populations and predicted reductions in favored food resources. Harris et al. (2006) projected how this population might fare in the future under alternative survival rates, and in doing so estimated the rate of population growth, 1983–2002. We address issues that remain from that earlier work: (1) the degree of uncertainty surrounding our estimates of the rate of population change (λ); (2) the effect of correlation among demographic parameters on these estimates; and (3) how a future monitoring system using counts of females accompanied by cubs might usefully differentiate between short-term, expected, and inconsequential fluctuations versus a true change in system state. We used Monte Carlo re-sampling of beta distributions derived from the demographic parameters used by Harris et al. (2006) to derive distributions of λ during 1983–2002 given our sampling uncertainty. Approximate 95% confidence intervals were 0.972–1.096 (assuming females with unresolved fates died) and 1.008–1.115 (with unresolved females censored at last contact). We used well-supported models of Haroldson et al. (2006) and Schwartz et al. (2006a,b,c) to assess the strength of correlations among demographic processes and the effect of omitting them in projection models. Incorporating correlations among demographic parameters yielded point estimates of λ that were nearly identical to those from the earlier model that omitted correlations, but yielded wider confidence intervals surrounding λ. Finally, we suggest that fitting linear and quadratic curves to the trend suggested by the estimated number of females with cubs in the ecosystem, and using AICc model weights to infer population sizes and λ provides an objective means to monitoring approximate population trajectories in addition to demographic analysis.
The grizzly bear (Ursus arctos) population inhabiting the Greater Yellowstone Ecosystem (GYE) has long been the object of both popular wonder and scientific scrutiny. Recently, we (Schwartz et al. 2006a) analyzed and summarized demographic characteristics of this population for 1983–2002. In generating projections useful for assessing effects of altered survival rates, we estimated mean λ at 1.042 or 1.076, depending on the treatment of radiocollared females whose fate remained unresolved (Harris et al. 2006). Despite treating a number of uncertainties in our analysis, some questions remained after publication of this work.
First, because our objective was understanding how a grizzly bear population with characteristics similar to the 1983–2002 Yellowstone population would behave in the context of differing levels of mortality (a prospective analysis, sensu Caswell 2000), we focused on the influence of process variance on λ (Harris et al. 2006). Sampling variance of demographic parameters was reported by Schwartz et al. (2006b,c) and Haroldson et al. (2006), but to the extent possible, excluded from the models of Harris et al. (2006). Thus, we provided no estimate of sampling error around our estimate of λ during 1983–2001, allowing some readers to erroneously interpret 1.042–1.076 as a confidence interval.
Second, our stochastic simulations assumed no process covariance among demographic parameters (Harris et al. 2006). One might reasonably wonder if this reflected reality, and if the absence of correlation structure in our models influenced our estimates of λ or estimates of sustainable mortality in the future.
Third, we showed that demographic data for a population roughly the size of the grizzly bears in Yellowstone — even were they somehow to be obtained in the absence of sampling error — can always yield misleading indications of the true underlying dynamic if analyzed over a short-time interval (3–5 years). Populations undergoing decline can produce occasional, short-term increases just as increasing populations can occasionally show short-term declines (Harris et al. 2006). We recommended analyzing demographics at roughly 8–10-year intervals (Schwartz et al. 2006d), but such a time span may be uncomfortably long for managers or the interested public. We also recommended continuing to index population size by estimating numbers of adult females in the population through mark–resight estimates of unduplicated females accompanied by cubs of the year. However, we did not make specific suggestions for interpreting trends using counts of females with cubs (given their inevitable stochastic fluctuations) because that was beyond the scope of our paper.
Here, we address issues raised above with additional analyses. First, in contrast to the work of Harris et al. (2006), we incorporate estimates of sampling variation to generate confidence limits on our earlier estimates of λ. Second, we use temporal covariates that were strongly supported in earlier models to indirectly estimate effects of correlation among demographic parameters on estimates of λ. Finally, we investigate methods to smooth annual fluctuations in estimated numbers of females with cubs in the population (Chao 1989, Keating et al. 2002, Cherry et al. 2007). We propose that information–theoretic methods be used to distinguish the current, linear increase from hypothetical future trajectories and estimate the power that such a monitoring protocol would have to detect these changes in population trajectory.
Study Area
The GYE, which we defined as our study area, encompasses all of Yellowstone and Grand Teton National Parks as well as portions of 6 adjacent national forests and smaller amounts of state and private land in Wyoming, Montana, and Idaho, USA. As of 2004, grizzly bears were estimated to occupy roughly 37,258 km2 (Schwartz et al. 2006e) of grasslands, shrub–steppe, open stands of juniper (Juniperus scopulorum), limber pine (Pinus flexilis), and Douglas-fir (Pseudotsuga menziesii), dense stands of lodgepole pine (Pinus contorta), subalpine stands of Engelmann spruce (Picea engelmannii), and whitebark pine (Pinus albicaulis), and alpine tundra and fellfields devoid of trees. Grizzly bears occurred at elevations of 1,584 to 3,656 m (Schwartz et al. 2002). For detailed descriptions of the study area see Blanchard and Knight (1991), Mattson et al. (1991), and Schwartz et al. (2002).
Methods
Field Methods and Derivation of Demographic Parameters
Demographic parameters were calculated from female grizzly bears radiocollared as part of research or management activities in the GYE, 1983–2001 (Schwartz et al. 2006a:11 discusses the two samples and how we obtained a representative sample). To estimate survival of independent females, 111 adult (age >5 yr) and 60 subadult (age 2–4) female grizzly bears were monitored via weekly or bi-weekly flights for various lengths of time (total 3,420 bear-months of monitoring, Haroldson et al. 2006). Survival was estimated on the logit scale using the known-fate (i.e., Kaplan-Meier) procedure in program MARK (White and Burnham 1999). Because some bears had unknown status at the time of last contact, 2 estimates were produced; one in which unexplained or unresolved losses were censored from the data at last contact (censored, C), and one in which all such animals were assumed to have died (all dead, AD). Cub and yearling survival rates were estimated by following the fates of 137 dependent young in 65 litters of 49 unique marked females. Schwartz et al. (2006c) used the nest success estimator in MARK to estimate survival on a daily basis on the logit scale. Reproductive rates were estimated from following 108 adult females (Schwartz et al. 2006b).
Reproduction and survival were modeled as functions of temporal and individual covariates (Schwartz et al. 2006a). Here, we used temporal covariates appearing in the best-supported models: an index to whitebark pine nut abundance (WBP); an index of winter severity (WSI), and an index of population size (Minpop). See Schwartz et al. (2006a,b,c) and Haroldson et al. (2006) for field methods, demographic analyses, and derivation of covariates.
Confidence Intervals Surrounding Estimates of λ
We used parametric bootstrapping (Alvarez-Buylla and Slatkin 1993, 1994) to calculate confidence intervals for λ. To model the 4 parameters Harris et al. (2006) used to estimate λ, i.e., reproductive rate (mx), cub survival (s0), yearling survival (s1), and independent (age >2) survival (sx), we developed beta distributions using parameters α and β, where
and μ was the mean and σ2 the variance (White 2000:298) of each distribution (e.g., sx, mx), as reported by Schwartz et al. (2006b,c and Haroldson et al. 2006; Table 1). We used the PopTools (G.M. Hood, 2004, PopTools version 2.6.2, http://www.cse.csiro.au/poptools) extension in Excel to run Monte Carlo iterations, sampling from all 4 distributions simultaneously, retaining the estimate of λ derived from that particular (randomized) combination of vital rates. For each iteration, λ was estimated using the Lotka formula where
, 29 was the maximum age considered, and λ = er (Caughley 1977). We ran 10,000 iterations for each of the 2 mean adult survival rates (AD = 0.922; C = 0.950; Haroldson et al. 2006), and interpreted the values between the 2.5th and 97.5th percentiles as the 95% confidence interval.
Table 1
Demographic parameters from Haroldson et al. (2006:35) and Schwartz et al. (2006b:20, c: 27), top row of each set, and realized rates from Monte Carlo simulations, bottom (n = 10,000) for a population of female grizzly bears with cubs in the Greater Yellowstone Ecosystem. Survival rates (sx) were estimated censoring (C) all unexplained or unresolved losses at last contact (censored, C), or assuming that all such animals had died (AD).
Correlation Structure of Demographic Parameters and Its Effect on λ
We used the top models of Haroldson et al. (2006) and Schwartz et al. (2006b,c) that employed temporal covariates to predict yearly survival of cubs, yearlings, and independent females, and of reproduction (summarized by mx). For mx, we used the model with the second lowest Akaike's Information Criteria (AICc) (ΔAICc = 0.58, Schwartz et al. 2006b:21, model 2), which included WBP and Minpop as temporal covariates (the top model did not include WBP). This model predicted litter size probability. We assumed constant inter-birth intervals and ages at first reproduction (neither of which can be estimated annually), and adjusted results until they produced mx of 0.318 (used in all previous projections) when both WBP and MinPop were at their mean values. For cub and yearling survival, we used the second highest ranking model (ΔAICc = 0.62; model 2, Schwartz et al. 2006c:29), which included only WSI. As an additional investigation, we used the third highest ranking model (ΔAICc = 0.85; model 3, Schwartz et al. 2006c:29) which included only MinPop. For independent female survival, we used the second highest ranking model (ΔAICc = 0.50), which included WBP and WSI (using the C data set; Haroldson et al. 2006 did not provide AD models because they were unduly influenced by the occasional presence of unresolved-fate bears). In each case, we used the yearly covariate (Schwartz et al. 2006a:14) and calculated the model's best estimate of each demographic rate for that year (with appropriate back-transformation and conversion to an annual rate). For mx (via litter size), we used the WBP and MinPop for year t − 1 to estimate the rate in year t (Schwartz et al. 2006b:13).
For both of these models, we ran 2 sets of Monte Carlo re-samplings of the 19 years of vital rates: (1) in which all rates varied in tandem, linked via their mutual association with that year's WBP, WSI, and MinPop values (i.e., co-varying rates using covariate values in Schwartz et al. 2006a:14), and (2) in which all rates varied independently, i.e., sampled from the 19 possible values, but scrambled from each other, in each case calculating asymptotic λ using PopTools and Eq. 3. A comparison of these 2 distributions (using identical vital rates, but one in which they were temporally correlated as suggested by our strongly supported models, and one in which rates varied independently) provided an estimate of the effect of co-varying rates on λ.
Interpreting Trends Using Estimates of Unique Females with Cubs
Under monitoring protocols (Interagency Grizzly Bear Study Team 2006) current as of 2006, an estimate is made annually of the number of females with cubs of the year (
), which are identified from all observations of females with cubs using a rule set developed by Knight et al. (1995). We consider sightings of females with cubs made with the aid of radiotelemetry biased and thus exclude them from the calculations (Keating et al. 2002). The Chao2 estimator (Wilson and Collins 1992, Keating et al. 2002)
f1
= number of bears seen exactly once,
f2
= number of bears seen exactly twice, and
m
= total number of unique bears seen,
is then used to estimate the total number of females with cubs present from the estimated number observed. This estimator, one of a number developed to deal with individual heterogeneity in capture probability, has lower mean squared error and percent relative bias than alternatives given the sampling intensity and recapture patterns observed in Yellowstone (Cherry et al. 2007). The trend in this segment of the population and its rate of change can be estimated from these annual estimates, providing an independent estimate of λ. However, annual estimates of 
can vary because of sampling error as well as because of synchronized reproduction. Thus, using annual estimates independently each year can result in greater variation in the estimate of total population size than is likely to characterize the true population.
Monitoring Population Size and λ Using Females with Cubs
The natural logarithm of the number of females with cubs [
] can be fit with a linear model of year (yi = i) to estimate λ as:
and the rate of population change is estimated as 
, giving 
. Asymmetric confidence intervals on λ can be estimated as the exponential of the symmetric confidence bounds on 
. Standard errors for ln( 
) can be computed with the usual linear model methods, and confidence intervals for (
) can be estimated as the exponential of the confidence bounds on ln( 
). When we assume a reasonably stable age and sex structure for the total population, this estimate of λ represents the rate of change of the entire population. Fitting a linear relationship makes the standard assumptions of least squares regression.
A quadratic regression can be also used to detect a change in 
through time. We fit the model
will become negative as population growth slows (as it would, for example, if the population reached carrying capacity) or reverses. Information–theoretic model selection methods (Burnham and Anderson 2002) can be used to select between the linear and quadratic models, and hence to detect changes in 
as additional data are collected. We used model averaging of the linear and quadratic models of the predicted population sizes of females with cubs and estimated λ to estimate population sizes through time (i.e., 
), and thus smooth the variation of the Chao2 estimates.
Power Analysis of Using 
to Estimate λ
To estimate the power of these data to detect a true reduction in λ (i.e., correctly choose the quadratic model), we estimated variance components of the Chao2 female-cubs counts 1983–2006, and applied these in Monte Carlo projections for 10 additional years under assumed values of λ.
To separate the sampling variance associated with each population estimate, 
from process variance, we fitted the linear model (above), with the assumption that the variance of 
was the sum of the sampling variance and the process variance. We found no evidence for serial correlation when regressing Chao2 estimates with time (Durbin-Watson D = 1.940). For the Chao2 estimator, 
was estimated with bootstrap resampling of the data, and the variance of the resampling distribution was the estimate of 
. We estimated the variance of 
using the Delta method as 
.
To estimate the process standard deviation from the 1983–2006 Chao2 estimates, we used PROC NLMIXED in SAS (SAS Institute, Inc. 2003). This procedure maximizes the likelihood of 
, and the process SD, with the likelihood specified as a normal distribution with mean predicted by 
and variance 
. This model thus explicitly includes the sampling variance of 
plus the process variance that is estimated by the procedure. Process standard deviation of 
was estimated to be 0.176 with SE = 0.0461 and 95% confidence interval of 0.0808–0.271
To estimate the expected sampling variance of future Chao2 estimates (which assumes that future sampling effort will remain approximately the same as that used to collect 1983–2006 data), the mean of the sampling variances of the natural log of the population estimates was computed. The expected sampling variance of future Chao2 estimates was then computed as a normally distributed random variable with mean of zero and standard deviation equal to the square root of mean sampling variance. From this procedure, the estimated sampling standard deviation was 0.34.
To evaluate sensitivity of the linear and quadratic models to changes in 
over 1 to 10-year intervals, we projected the 2006 population estimate of N2006 = 52 (obtained by model averaging the linear and quadratic model estimates from the fits of the 1983–2006 data), assuming alternative λ values of 0.95, 0.975, and 1.0, and using our estimates of process and sampling variation (above). Population size for each succeeding year was generated with the recursive relation 
, where the process variation was added as 
, a normally distributed random variable with mean zero and standard deviation 0.176. The estimated population size (corresponding to the Chao2 estimates) was taken as 
, where the sampling variation 
was added as a normally distributed random variable with mean zero and standard deviation 0.34. Each replicate was simulated independently: completely new data were added to the 1983–2006 data for each simulation.
One thousand replicates of each of the 30 scenarios (3 alternative λs × 10 alternative time-frames) were generated, from which we estimated the mean AICc weight of the quadratic model, the proportion of iterations in which the quadratic term was selected (weight > 0.5), and the power of the t-test to reject the null hypothesis that the quadratic term was equal to or greater than zero. This realistically simulated the data and analyses that would be available to managers when judging whether the population had changed its trajectory in a downward direction.
Results
Confidence Intervals Surrounding Estimates of λ
Survival and reproductive rates reported by Haroldson et al. (2006) and Schwartz et al. (2006b,c) were similar to those realized from our Monte Carlo re-sampling (Table 1). Under the AD assumption, the approximate 95% confidence interval around the 1983–2001 λ of 1.041 was 0.972–1.096, and about 10.3% of re-samplings suggested a negative trend. With all such animals censored at last contact (assumption C), the approximate 95% confidence interval around the estimated λ of 1.066 was 1.008–1.115, and less than 2% of re-samplings suggested population decline.
Correlation Structure of Demographic Parameters and Its Effect on λ
As expected, survival of adults, cubs, and yearlings were positively correlated (Table 2). (The correlation between yearly cub and yearling survival was imposed by the model of Schwartz et al. 2006b). In contrast, correlations among reproduction and survival were weak and inconsistent (Table 3). Although WBP appeared in all models (and the direction of the relationship was always the same), the reproduction model used WBP in the previous year, whereas survival models used WBP in the year of survival.
Table 2
Point estimates of yearly survival and reproduction of Greater Yellowstone Ecosystem grizzly bears as predicted by the second ranking models; i.e., those using temporal covariates; mx values lag one year behind. WSI is an index of winter severity; Minpop is an index of population size.
Table 3
Correlation (above diagonal) and covariance (below diagonal) matrix among yearly rates predicted by models described in text for a population of female grizzly bears with cubs in the Greater Yellowstone Ecosystem.
When using cub and yearling survival rates from the WSI-only model, Monte Carlo re-sampling (n = 5,000) with the implied correlation structure yielded a mean λ of 1.0652 (95% CI = 1.0084–1.1276). Identical re-sampling with rates scrambled independently of one another (n = 5,000) yielded a mean λ of 1.0648 (95% CI = 1.0182–1.1257). When using cub and yearling survival rate estimates from the MinPop only model, Monte Carlo re-sampling (n = 5,000) with the implied correlation structure yielded a mean λ of 1.0720 (95% CI = 1.0063–1.1335). Identical re-sampling with rates scrambled independently of one another (n = 5,000) yielded a mean λ of 1.0713 (95% CI = 1.0245–1.1314).
Interpreting Trends Using Estimates of Unique Females with Cubs
Monitoring Population Size and λ Using Females with Cubs
We used Chao2 estimates for 1983–2006 (Table 4) to estimate the rate of population change (Fig. 1). The parameter estimates and AICc weights for the linear and quadratic models (Table 5) suggested that primarily the linear model was needed to model changes in the number of females with cubs during this period. The estimate of λ using the linear model was 1.046 with 95% confidence interval of 1.030–1.061. The estimated quadratic effect (–0.000968, SE = 0.0012) was not significant (P = 0.419), and 74% of the AICc weight was associated with the linear model. Thus, the linear model was the best approximating model for 1983–2006, but we also provide the model averaged estimates (Fig. 1).
Table 4
Population estimates of female grizzly bears with cubs from the Chao2 estimator for 1983–2006, the Greater Yellowstone Ecosystem.
Figure 1
Model-averaged estimates of the number of female grizzly bears with cubs 
for 1983–2006 (connected points), where the linear and quadratic models of 
were fitted. The solid line represents the model-averaged slope of 
on time. The inner set of dashed lines represents a 95% confidence interval on the predicted population size, and the outer set of dashed lines represents a 95% confidence interval for the individual population estimates.
Table 5
Estimates and model selection results from fitting the estimates of female grizzly bears with cubs from the Chao model for 1983–2006, Greater Yellowstone Ecosystem.
Power Analysis of Using 
to Estimate a Change In λ
When our best estimates of process and sampling variation were added to hypothetical years 1 through 10, approximately 5 years were required of the population decreasing 5% yearly (i.e., λ = 0.95) before the preponderance of evidence (AICc weight > 0.5) favored the quadratic model (i.e., fundamental change in state from linear increase, Fig. 2). Under the scenario in which population size stabilized after year 2006 (i.e., λ = 1.0), 7 or 8 years were required for the preponderance of evidence to favor the quadratic model (depending on the criterion used, Fig. 3). Power to detect a yearly decline of 2.5% was intermediate between these 2 examples.
Figure 2
Mean AICc weight of the (negative) quadratic term, proportion of simulations in which the quadratic model had greater AICc weight than the linear model, and power of the quadratic term (i.e., probability of rejecting the linear model) when expected λ changed to 0.95 following the 1983–2006 series of estimates of females with cubs, for additional years 1 to 10 and using estimates of process and sampling variation from the data.
Figure 3
Mean AICc weight of the (negative) quadratic term, proportion of simulations in which the quadratic model had greater AICc weight than the linear model, and power of the quadratic term (i.e., probability of rejecting the linear model) when expected λ changed to 1.0 following the 1983–2006 series of estimates of females with cubs, for additional years 1 to 10 and using estimates of process and sampling variation from the data.
Discussion
Confidence Intervals Surrounding Estimates of λ
We observed general concordance in our estimates of λ during 1983–2001 using demographic data (x̄ = 1.041, 95% CI = 0.972–1.096 assuming all radioed females with unresolved fates died; x̄ = 1.066, 95% CI = 1.008–1.115 with those females censored at last contact) and of λ during 1983–2006 using the Chao2 estimator applied to counts of females with cubs (x̄ = 1.046, 95% CI = 1.030–1.061). Because the raw data underlying these 2 approaches are independent, we believe this concordance provides confidence that the apparent increase in the population was real. That said, because each approach makes different assumptions (and in this case, cover differing time spans), a strict comparison is impossible. In the case of the demographic model, our inability to resolve the fates in the field of a few adult females produced uncertainty equal to approximately a 2.5% yearly rate of change, and under either assumption, confidence intervals were broad. In the case of counts of females with cubs, even the most rigorous of field and office protocols leave room for uncertainty about m, f1, and f2, and in general, the method of identifying unique females with cubs (Knight et al. 1995) inevitably yields an increasingly conservative estimate of m as the density of adult females increases.
Our Monte Carlo resampling for obtaining a confidence interval around λ is similar to the bootstrapping method used by Eberhardt et al. (1994) in their analysis of Yellowstone grizzlies (Hovey and McLellan 1996, McLoughlin et al. 2003, Garshelis et al. 2005). We used beta distributions (Powell et al. 2000) rather than resampling the raw data because we had large sample sizes. Given the influence of a few females whose fate remained unresolved, we believe small differences in the method used to obtain confidence intervals around λ are inconsequential and that our nonparametric bootstrapping provides adequate representation of the confidence interval (Alvarez-Buylla and Slatkin 1994).
Correlation Structure of Demographic Parameters and Its Effect on λ
Ignoring yearly correlation among rates had little effect on estimates of the mean λ (e.g., λ = 1.0648 assuming complete independence, λ = 1.0652 with positive correlation among parameters). However, the width of 95% confidence intervals surrounding λ generated in the absence of covariance in demographic parameters may be slightly underestimated. We lack a method to quantify the practical effect of this underestimation directly; our approach here explicitly included sampling variation, whereas the simulations of Harris et al. (2006) reflected shrinkage estimates of survival for independent females because including sampling variation was inappropriate in a projection model. Thus, we were not surprised to find wider confidence intervals. Given that the best supported models of Haroldson et al. (2006) and Schwartz et al. (2006b,c) suggested some positive correlation among demographic rates, we recommend that the confidence intervals (assuming independence among parameters) we report above (0.9723–1.0961 and 1.0076–1.1154), as well as the estimates of variability around λ given specified mean survival and levels of temporal variation (i.e., Harris et al. 2006:50; Table 20) be viewed as being slightly narrow.
Interpreting Trends Obtained Through Monitoring Unique Females with Cubs
Adult females are the critical segment of the population: they control reproduction, and their survival rates are the most important single component of overall population trajectory (Eberhardt et al. 1994, Hovey and McLellan 1996, Garshelis et al. 2005, Harris et al. 2006). Thus, estimating the rate of change of this segment as one measure of the rate of change of the entire population is appropriate. The Chao2 estimator remains our choice as the best way to estimate the total number of females with cubs yearly from those observed (Cherry et al. 2007).
Unfortunately, the number of females with cubs can only be estimated, and even our best field and analytic procedures can only do so much to reduce variation in yearly estimates (although any reduction in sampling variance would increase power to detect trends). We do not doubt that these estimates track true population change generally, but they have limited power to detect subtle (yet potentially important) changes in λ within a few years.
We suggest fitting linear and quadratic models to the time series and using AICc weights to evaluate their relative appropriateness rather than using a traditional hypothesis-testing framework. Although statistically sound, significance tests implicitly place the burden of proof on the hypothesis that change has occurred (the quadratic model in this case), which, in this situation, could easily lead to a long delay before data are sufficient to reject the null. Further, hypothesis testing results in an either–or situation, where inferences are made from only 1 of the 2 models. In contrast, model weighting under the information–theoretic paradigm weights the 2 models relative to their respective likelihoods. That said, we note that even using the AICc weight criterion given realistic projections under variability, a delay of some years is inevitable before one gains appreciable power to correctly detect a change of the magnitude we modeled. This degree of uncertainty reflects not so much a poor choice of index or analytical technique as it does the inevitable characteristics of field work on grizzly bears and the mathematically (if not biologically) subtle distinction between growth at about 5% annually and a stable or slightly declining population.
Management Implications
Managers and the interested public would like to know not only the growth rate of Yellowstone grizzly bears in the past, but how the population is changing now. Alas, we believe that the quest for absolute certainty or a quick answer is likely to be futile. While we feel confident that the GYE grizzly bear population increased substantially during 1983–2006, we acknowledge that the confidence interval around estimates of λ derived from counts of females with cubs is broad, and that from our demographic analyses, even overlaps 1.0 (under the most conservative assumptions). Nonetheless, the close match between the 2 analyses from independently obtained data sets (as well as corroborating evidence such as the evident increase in range, Schwartz et al. 2006e) increases our confidence that the population has grown.
We recommend continued monitoring of females with cubs of the year, fitting both linear and quadratic models to the complete data set, and using AICc to evaluate the strength of these 2 competing models. Weight favoring the quadratic term is evidence that population growth has slowed or reversed, but we caution that lack of such evidence is not necessarily proof that change has not taken place. Gradually increasing evidence for the quadratic model over a few years (assuming a negative quadratic slope) should serve to keep biologists and managers alert to a possible change in system state. If the AICc weight favors the quadratic term (i.e., >0.5) in modeling the rate of change of females with cubs, we recommend that a biology and management review (USFWS 2007) be undertaken. Under the best of circumstances, this monitoring protocol leaves uncertainty about the system state during the most recent few years. We find this compelling reason to couple counts of females with cubs to continued monitoring of demographic rates from a sample of radiomarked females and their offspring. Although also characterized by variability and time-lags, such monitoring provides an independent measure of population vigor and is likely to be helpful in explaining any observed changes in numbers of females with cubs.
Acknowledgments
We thank S. Cherry, K. Keating, D. Moody, and C. Servheen for discussion and review. Support for this work was provided by US Geological Survey Biological Resource Discipline. We thank R.R. Knight and B.M. Blanchard for their efforts in data collection and program direction with the study team from its inception in 1973 through 1997. We also thank the many observers over the years who faithfully provided us with sightings of females with cubs of the year in the Greater Yellowstone Ecosystem. Funding for counts of females with cubs was contributed by members of the Interagency Grizzly Bear Study Team, including the US Geological Survey, Northern Rocky Mountain Science Center, Wyoming Game and Fish, Montana Fish Wildlife and Parks, Idaho Fish and Game, the National Park Service, US Forest Service, and the US Fish and Wildlife Service. Without the continued support from all agencies, this monitoring program would not have been possible.
